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The original goal of Large Eddy Simulations of fully developed turbulent flows was to accurately 

describe large-scale flow features u(A) at the scales r > A whore A is a size of computational mesh. 
The effect of small-scale velocity fluctuations (r < A) was to be accounted for by effective transport 
coefficients (subgrid models) in the coarse-grained Navicr-Stokcs equations. It is shown in this paper 
that, due to anomalous inertial range scaling (intermittency) of the moments of velocity difference, 
the existing subgrid models are intrinsically incapable of quantitatively describing flow features at 
the scales r < NA with N w 10. This increases computational work approximately by a factor 
10^ — 10*. The breakdown of the widely used Smagorinsky relation for the subgrid viscosity on the 
scales A/L < 1 is demonstrated and a modification accounting for intermittency of the filtered out 
small-scale fluctuations is proposed. 



Introduction. In 1969 S. Orszag introduced spectral 
methods as a tool for direct numerical simulations (DNS) 
of strong turbulence [1] which were implemented in simu- 
lations of three - dimensional isotropic flow in 1970-1971 
and published in his paper with Patterson in 1972 [2]. 
Thus, a new field of numerical experimentation on tur- 
bulent flows was born. Based on Kolmogorov's theory 
Orszag estimated the scaling of computational work W 
with the large-scale Reynolds number Re = UrmsL/vo'- 
W = 0{Re^). Here L and uq stand for the integral scale 
corresponding to the top of inertial range and molecular 
viscosity, respectively. Although it became clear recently 
[3] that, due to intermittency, the scaling of computa- 
tional work required for full simulations of developed tur- 
bulence can be as bad as W = 0{Re'^), this work, defln- 
ing computation as an alternative experimental tool for 
turbulence studies, was a major breakthrough. 
It became clear immediately that due to the unfavorable 
Reynolds number scaling of computational work, the 
DNS of high Reynolds number flows of engineering in- 
terest were impractical and derivation of coarse-grained 
equations (models) for prediction of the large-scale flow 
features, theoretically accounting for the small-scale 
fluctuations were of crucial importance. The aim was to 
eliminate the small-scale velocity fluctuations, describing 
structures on the scales r < A. from the Navier-Stokcs 
equations, express their effect in terms of the resolved 
field and reduce the computational work toW = 0{Re'^) 
with 7 w — 1. 

The problem is formulated as follows. The Reynolds 
stress can be decomposed as: 

RS = u • Vu + u • VusG + usG • Vu -h Msg ■ Vusg 

where u and Msg a-re the resolved and subgrid velocity 
fields describing flow structures on the scales r > A 



and r < A, respectively. To fully simulate velocity 
fluctuations on the length-scale r varying in the interval 
L > r > riK ^ LRe~ i , one needs a computational mesh 
at least as small as A « ry^- If due to computational 
limitations the practically available mesh size A ^ rjn, 
one needs equations for velocity fleld u(A) = u defined 
on the resolved scales r > A. The rapidly varying terms 
involving the computationally unaccessiblc small-scale 
subgrid fluctuations Usg on the scales rjx < i" < A are 
to be averaged out and their effect represented in terms 
of the resolved field u. 

First attempts to achieve this goal were made by Dear- 
dorff [4] who applied the model previously derived by 
Lilly and Smagorinsky [5] to simulations of wall flows. 
This model is based on an assumption that the effect 
of small - scale velocity fluctuations on the large-scale 
ones can be represented in terms of effective viscosity. It 
will become clear below that while this assumption is ac- 
curaie when A/L w 1, due to the anomalous scaling of 
the inertial range structure functions, it is grossly incor- 
rect when the cut-off A is small, i.e. A/L ^1. The 
essence of the Lilly-Smagoronsky construction can be for- 
mulated as follows. Consider a length - scale A in the 
inertial range. Then, according to Kolmogorov's theory, 
the scale-independent energy flux 

E = f (A) « 2v{A)S^{A) « 
\{u{x + A) - u{x)Y\/ A = ^3(A)A2 (1) 

where 5(A) = y/Sij{A)Sij{A), with ^^^(A) = ^{diUj -h 
djUi), is the rate of strain evaluated on the resolved 
scales r > A and the effective viscosity z^(A) accoimts 
for the contribution of the small-scale velocity fluctu- 
ations acting on the scales r < A. The filtering of 
the small-scale velocity fluctuations from the interval of 
scales r < A generates the resolved field which is smooth 
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on the scales r > A and, as a result, the spatial deriva- 
tives djUi{A) ft! (Mi(x + Ajj) — u{x.))/Aj are expressed 
in terms of velocity increments of the scale A. Prom 
here, we have for the well-known relation for turbulent 
viscosity viA) ft A2|S'(A)3|/52(A). 
This global argument, dealing with the mean proper- 
tics of a flow, is relatively safe. To develop a numerical 
method, one needs local differential equations for the nu- 
merically resolved velocity field with a local expression 
for turbulent viscosity in terms of this resolved field. A 
dimensionally correct possibility is the following : based 
on (1), write 5(A) = 2u{A) Sij (A) Si j (A) + ^{A) where 
^' is an unknown function satisfying a global constraint 
\E' = 0. Evaluation of this function which is an extremely 
difficult task will be discussed below. The simplest and 
crudest way to develop a model is by dropping the av- 
eraging symbols in (1), setting \E' = and, interested in 
isotropic turbulence, taking for simplicity Aj = A obtain: 



the resolved velocity field (for state-of-the-art review, see 
Ref.[7]) 

dtu + WiV^u = -Vp + VaJ^s(A)VaU + EN + { (4) 

where z^(A) originates from the usg ■ Vu^g contribu- 
tions to the Reynolds stress and EN stands for the back 
scattering term ("eddy noise") leading to a small frac- 
tion of energy the large-scale eddies receive from inter- 
actions with the small-scale ones. Interested in three- 
dimensional turbulence dominated by the direct cascade 
we, for a time being, can neglect the EN- contribution. 
The large-scale forcing function f in the right side of (4) 
is not modifled by the coarse-graining (filtering) proce- 
dure described below. Multiplying (4) by u, the energy 
balance for the resolved velocity field is calculated read- 
ily: 



us{A) =aSA'^ K a\u{x + A)-u{x)\A (2) 

which is basically a mixing length model well-known 

from Prandtl's works. The derivation of this model di- 
rectly from the Navier-Stokes equation, based on the one- 
loop renormalized perturbation expansion led to the ex- 

_ A 

pression for effective viscosity Veff = J^s(l — e ''^f ) 
where t]k ~ LRe~i is the Kolmogorov dissipation scale 
[6]. Typically, in the LES the cut-off A » rjK, so that 
i^ef f ~ i^s- This expression for effective viscosity, com- 
bined with the Navier-Stokes equations for the resolved 
velocity field, is a celebrated Smagorynsky model for the 
large-eddy-simulations of turbulence (LES). The expres- 
sion (2) is the basis of majority of modern LES schemes 
postulating the length-scale A identical to the size of 
a cell of computational mesh. We can see that since 
i/(A) >> Vq, the role of this increased viscosity is to re- 
duce the Reynolds number and, consequently, the num- 
ber of degrees of freedom necessary for description of the 
large-scale flow features. In general fluctuating param- 
eter a, which is constant in isotropic and homogeneous 
turbulence cannot be derived from this qualitative argu- 
ment. 

As follows from (2) the moments of fluctuating resolved 
dissipation rate are: 

£^^={usS^)^ oc a"|(^A«|)3")/A" oc (£L)S (^)«^"-" 

(3) 

with "anomalous" scaling exponents ^„ consistent with 
the so called Kolmogorov's refined similarity hypothesis. 
The relation (2) is a result of a general small-scale 
filtering procedure or, equivalently, the scale elimina- 
tion methods developed for the dynamic renormalization 
group. All existing sub-grid models are based on an as- 
sumption that the filtering leads to the equation for 



V = ^f = aA'\S\{^y = 

^. \u{x + A) - u{x)\{u{x + A) - uix))"^ . 
^ A ' 

which, by choosing a proper magnitude for the constant 
a, is consistent with the Kolmogorov's 4/5 law. Since 
by construction, the arbitrary cut-off A is chosen in the 
inertial range, the equation (4) with the Smagorinsky vis- 
cosity (2) is expected to produce approximately accurate 
results for the resolved energy spectrum and the second 
moment S2 = {u{x + A) — u(x))^. 

However, if consistent, the local equation (4) with sub- 
grid model (2) is supposed to describe the resolved ve- 
locity field, including all structure functions Sn{r) = 
{u{x + r) — u{x))^ and moments of resolved dissipation 
rate on all scales r > A. 

Scale by scale comparison of LES with DNS . To assess 
the scale-by-scale performance of a subgrid model we 
have introduced the displacement r = A''A where N > 1 
are the integers, and evaluated the modified structure 

functions 6'p(A''A) = \uix + N A) - u{x)\p'' = \Snau\p'' 
using the DNS (resolution 256^ ) and variational LES 
(resolution 32'^) described in detail in the Ref. [8]. For 
p = 2 this parameter gives ^/S2 = (^jva w)rm« and for 
an arbitrary p it is a useful tool for analyzing quality of 
prediction of the most intense small-scale events. The 
energy spectra E{k), evaluated using both methods in 
Ref. [8], were in a good agreement with each other [8]. 
On Fig.l we show the ratio Sp^LES I Sp^oNS as a function 
of displacement r = NA. As expected, even for this low 
Reynolds number flow for which A/L ft 0.1, the "small- 
scale" results rapidly deteriorate with increase of power 
p. This drawback of all existing subgrid models (for dif- 
ferent examples see an excellent review [7] ) , which can- 
not be cured by introducing the scale-dependent 
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FIG. 1: Testing subgrid model: = '''"'^"'"^i' 
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algebraic relations a = a(r/A), is quite costly: the 
model which is capable of accurately describing the flow 
on the scales from the interval r/A > 10 only, is at least 
by the factor 1000 more expensive than the one valid in 
the entire interval L > r > A. 

The Smagorinsky model has been derived in Ref. [6] 
directly from the Navier-Stokes equations from the one- 
loop renormalized perturbation expansion in powers of 
the Reynolds number Re^ff = u^msL / viL) = 0(1) 
which can be reformulated as a series in powers of di- 
mensionless rate -of-strain proportional to the velocity 
increments with displacement r in the inertial range. As 
well-known from experimental data, at the large scales 
r ^ L, the velocity fluctuations obey Gaussian statistics 
and, therefore, all contributions to the expansion are of 
the same order and a few flrst contributions to the series 
give quantitatively accurate results. As r/L << 1 this 
is not so high-order-contributions to the series become 
dominant and must be taken into account. The sever- 
ity of this problem, stressed in [9] , becomes clear from a 
general expression for the probability density of velocity 
fluctuations (increments) [10]: 

(5) 

where, to simplify notation, we set r = r/L and SrU = 
Sru/urms- In the limit r — >■ 1, this expression gives Gaus- 
sian PDF P (X /2 — which is a probabil- 
ity density of single-point velocity field. 
The small-scale {r/L < 1) strong luctuations are well 
described by the expression following from the steepest 



descent evaluation of the integral (5) (for details see [10]- 
[11]): 

2 r°° 2 (ln^^)2 
p(S u) = - / e^"^ Exn\~- '"^^^ ]d- 

(6) 

The constants in (5)-(6) tested in various flows [11]-[13], 
are a « 0.383 and b « 0.01667. As A/L ^ 0, the 
PDF (6) is an extremely shallow function stressing 
the importance of strong small-scale fluctuations not 
accounted for by Smagorinsky-like models. Thus, to 
eliminate small-scale fluctuations r < A, one has to use 
(5) -(6) as a fllter applied to the Navier-Stokes equations 
which is not a simple task. 

To illustrate the role of intermittency (anomalous scal- 
ing) in the LES modeling we, can use an approach 
leading to derivation of turbulence models from a ki- 
netic equation for some artificial "particles" forming a 
"fluid" characterized by the relaxation time r reflecting 
dynamic properties of a system. For example, choosing 
T = const leads in the first order of the Chapman-Enskog 
expansion to the Navier-Stokes equations with viscosity 
i^Q = Tu^^g/i with non-linear corrections appearing from 
the higher order terms. The same procedure but with the 
relaxation time t oc JC/£ yields the well-known linear and 
non-linear K. — £ models for the large-scale velocity field 
[14]. Repeating the procedure, presented in Refs.[14], it 
is easy to see that the LES, Smagorinsky-like, models are 
generated by the same equations with t(A) = l/\Sij{A)\. 
In the second order of perturbation expansion in powers 
of dimensional rate of strain the procedure gives: 

a,, = 2iys{A)S,, - 2iys{A){dt + u • V)(t5,,)/3 
-4:1^3 {A)T[SikSkj - 2^ijSkiSki] + 

iysiA)T[Sik^k] + Sjk^kj] (7) 

with i^s = Tu2,„^(A) = = aA'^\S\. The same re- 

lation has been derived by Rubinstein et al [15] directly 
from the Navier-Stokes equations using scale-elimination 
technique and, even earlier, by Speziale [16] who ap- 
plied symmetries used in rheology for derivation of tur- 
bulence models. On the first glance dropping the aver- 
aging symbol from expression Tu'^^g{A), which is used 
for derivation in the procedure leading to Smagorinsky 
model, cannot be accurate. However, a consistent deriva- 
tion leading to the expression (7) involves averaging over 
relatively fast velocity fluctuations varying on the time 
scales rir < A) < t(A). Thus the velocities on the scale 
r « A can approximately be considered "slow modes". 
This observation somewhat justifies an assumption lead- 
ing to (7). As we see, the Smagorinsky model, which 
is the basis of all modern numerical schemes for LES is 
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obtained in the lowest order of expansion in powers of 
dimensionless rate of strain. Now, we will discuss the 
limits of its validity. 

Our goal is to express the s\ib-grid contribution to the 
Reynolds stress is terms of the resolved field u(r > A). 
As one can see from (7), due to prolifiration of the sub- 
scripts, with increase of the power of expansion each 
term becomes more and more comlex. Some simplifi- 
cations are based on the following: The resulting equa- 
tion must be invariant under Galileo transformation. 
Only two dimensionless operators satisfy this require- 
ment: Di = T{dt + UiVi), responsible for the memory 
effects, and D2 = t^^. Thus, the general expression for 
the Reynolds stress is aij — aij{Di, D2). 
It has been shown [17]-[18] that expansion in powers of 
Di results in the telegrapher equation with 0( rZ)i)-term 
in the left side of (4, ) responsible for the quasielastic re- 
sponse of turbulent flow to external high-frequency per- 
turbations. In what follows, mainly interested in the in- 
termittency contributions to (Tij , we omit this term which 
is important for the flows rapidly varying on the time 
scale T << t(A). To simplify the treatment and to avoid 
dealing with an extremely complex tensorial structure of 
high-order terms in the expansion (7) we, based on the 
constraints discussed above, choose an infinite subset of 
leading contributions to each term in power series of aij 
similar to the " single-back-bone" diagrams of the Wyld's 
expansion of the Navier-Stokes equations: 

..,«.s5.,(A)^a„(§4§)" (8) 

where the coefficients Q!„ can be evaluated only for the 
first few terms. Similar approach has been developed for 
the scalar problem in Ref.[19], where it was used for as- 
sessment of the accuracy of the e-expansion in turbu- 
lence theory. It is clear that the mean dissipation rate 

£ = —Ui^-^^ and the sub-grid viscosity can be defined 
as: 

...(A)=,.,5-..,i:<.„(M||)" ,9) 

This expression is easily understood if we notice that the 
expansion we are dealing with is in powers of dimen- 
sionless rate-of-strain rj = t{A)S similar to that discov- 
ered in derivation of the RNG K. — £ model from the 
Navier-Stokes equations [20]. Adding to this expression 
contributions including fljj does not change qualitative 
conclusions of the argument. 

As we see, in case of normal scaling or in the large- 
scale limit AIL 1, the ratios ( ) = B„ are 
A-independent and the dissipation rate based on this 



sub-grid viscosity in the limit A ^ L is proportional 
to 



n 

provided OLnl^n is a finite number. This expression is 
equivalent to Smagorinsky model (2). In the large-scale 
limit A — > L — > 00, due to the Galileo invariance, the 
corrections to the non-linear term u • Vu in the coarse- 
grained equation (4) are equal to zero and, as A/L w 
1, the energy balance in the equation (4) with viscosity 
given by (9) is not violated. 

A completely different result is obtained when anoma- 
lous scaling is taken into consideration on the; scales 
A < L. It is clear from (9) that in this case the dissipa- 
tion is not equal to (10) but is multiplied by an infinite 
sum: 

£{A) « vsg{A)S^{A) « f (L) ^a„/3„(^)«"+3-t-i « V 

n 

(11) 

Since ^„+3 - | - 1 < 0, for A/L << 1 the dominat- 
ing contributions from the high-orde terms in (10) in- 
validate the Smagorinsky's model in the range of scales 
A/L « 1. As we see, since production P is a large- 
scale A-indepcndent property, when A/L << 1, the en- 
ergy balance in the LES equation (4) can be restored 
only by renormalization of the non-linear term, i.e. by 
accounting for the high-order non-linear contributions, 
disappearing in the limit A/L —> 1. In this case the rela- 
tion proportional to (9) is a possible candidate satisfying 
the above requirements. Thus, neglecting the "rapid- 
distortion" (memory) contribution, the LES models ac- 
counting for the inertial range intermittency is: 



dt\i + ^{S, -KViU = -Vp-hV«i^sG(A)VaU-|-£iV-hf 

(12) 

As A/L ^ 1 the functional $ ^ 1. In the interval 
A < L it ensures cancellations leading to the cut-off 
-independent expression for the mean dissipation rate 

£{A)=£ = V. 

To conduct numerical simulations using equations with 
transport coefficients given by infinite series in powers 

of dimensionless ratc-of-strain tS is not a simple task. 
One can attempt to qualitatively represent the sum in 
(9) in a compact form like vsg{A) = vs/{l + =2) which, 

by the averaging the expression for the dissipation rate 
using the probability density (5)- (6) is equivalent to the 
infinite series (11). This qualitative representation of an 
infinite series (9) proved to be very useful in the K, — £ 
modeling emerging from a present theory in the limit 
t( A) ^t(L) oc/C/f [21]. 
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To summarize: in this paper we compare the LES based 
on Smagorinsky sub-grid model with the direct simula- 
tion of the same flow. It has been demonstrated that, 
while the LES gives a reasonably accurate representa- 
tion of the large-scale flow features, it completely under- 
estimates strong small-scale (r < 10 < A) intermittent 
velocity fluctuations. The reason for this failure is traced 
to the inadequate account for the small-scale intermit- 
temcy of the filtered out small-scale fluctuations in the 
Smagorinsky-class models. An improved model including 
the intermittent phenomena is proposed. 
It is interesting to compare the outcomes of this work 
with the results presented in two important papers [22]- 
[23]. In these works the authors compared filtered and 
unfiltered experimental data on turbulent properties of 
a few flows (wakes, grid, decaying turbulence, etc). The 
flltering procedure of Refs.[22]-[23] was the same as the 
one used in the LES. Since no modeling was involved, and 
one would expect that the resolved velocity fluctuations 
on the scales r < A of the two fields must be identical. A 
completely different picture emerged: the strong small- 
scale features of the filtered field were grossly underpre- 
dicted compared with the raw, unfiltered, data. This im- 
portant result indicates that large and small -scale veloc- 
ity fluctuations strongly interact and one cannot obtain 
a reasonable representation of one by simply filtering out 
the other. This conclusion becomes clear if one recalls 
that in turbulent flows extremely thin and long geomet- 
rical structures, like pancakes, "worms" etc, consist of in- 
separable, strongly interacting large and small-scale flow 
features. Thus, smoothing the velocity field over length- 
scale A (filtering), distorts the entire field including the 
interval A < r < lOA. In short, the works [22]-[23] 
showed that the simple-minded filtering procedure is in- 
trinsically fiawed. The theory described in this paper, 
based on consistent procedure of small-scale elimination 
from the Navier-Stokes or Boltzmann-BGK equations, 
proposes a remedy for this drawback by strong modifi- 
cation of the equations of motion for the resolved scales. 
The resulting model, while at the large scale equivalent 
to the ordinary Smagorinsky model, at the small scales 
r i=a A does not resemble neither it nor the Navier-Stokes 
equations themselves . 

Correctly accounting for the entire range of velocity 
fluctuations L > A > tjk may substantially reduce cost 
of numerical simulations and improve quality of LES of 
mixing and sound generation by turbulence. This will 
be a subject of future communications. 
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